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Abstract 


This paper is focused on the Bidirectional Reflectance Distribution Function (BRDF) in the context of algorithms 
for computational production of realistic synthetic images. We provide a review of most relevant analytical BRDF 
models proposed in the literature which have been used for realistic rendering. We also show different approaches 
used for obtaining efficient models from acquired reflectance data, and the related function fitting techniques, 
suitable for using that data in efficient rendering algorithms. We consider the algorithms for computation of 
BRDF integrals, by using Monte-Carlo based numerical integration. In this context, we review known techniques 
to design efficient BRDF sampling schemes for both analytical and measured BRDF models. 


Categories and Subject Descriptors (according to ACM CCS): 1.3.7 [Computer Graphics]: Three-Dimensional 
Graphics and Realism. Color, shading, shadowing, and texture. I.6.8 [Types of simulation]: Monte Carlo 


1. Introduction 


One of the goals in Computer Graphics is the production of 
synthetic images from digital object models. For some ap- 
plications when photorealism is required, those images must 
exhibit some degree of visual realism. Here, the term visual 
realism means that human observers should notice as little 
difference as possible between a synthetic image and a real 
photograph of the real object whose digital model was used 
to produce the image. 


This realism can be achieved by the production of a de- 
tailed model of geometry and materials in the scene. Part of 
the effort devoted to modeling must be focused on modeling 
the reflective properties of materials, as this is a key factor 
in visual realism. This is done by characterizing the direc- 
tional distribution of reflectance of the materials using phy- 
sical models. These models must be suitable for the design 
of efficient global illumination computation algorithms. 


The characterization of reflective properties of materials 
can be achieved by defining the function which determines 
how reflected radiance is distributed in terms of the distribu- 
tion of incident radiance. This function is the Bidirectional 
Reflectance Distribution Function (BRDF). 


This paper is meant to serve as a tool for researchers and 
developers. We review the Computer Graphics literature re- 
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lated to the BRDF function, from theoretical models to prac- 
tical BRDF-related aspects of implementation of rendering 
software. 


The paper is organized as follows: in Section 2 we provide 
the formal definition of the BRDF, and we introduce its pro- 
perties and main types. In Section 3, we include a review of 
main analytical BRDF models proposed in the literature. For 
each model we show their properties, advantages and disad- 
vantages in terms of realism and computational efficiency. 
Section 4 is devoted to acquired BRDF models. In Section 
5, we focus on the methodologies for BRDF sampling in the 
context of Monte-Carlo based numerical schemes for radi- 
ance integration. Section 6 finalizes with the conclusions. 


2. Definition and properties of the BRDF function 


Generation of realistic images relies on the numerical com- 
putation of approximations to the radiance function (denoted 
by L). The importance of L function in rendering comes 
from the fact that this quantity can be used to create sim- 
ple and accurate models of image formation in the human 
visual system and in photographic cameras [Shi90]. In this 
section we introduce formal definitions and properties for 
several radiance-related functions. A valuable introduction 
to concepts and terminology can be found in [NRH*92]. A 
far more detailed introduction (from the perspective of trans- 
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port theory) can be found in [Gla94], while more practical 
approaches (specifically focused on Computer Graphics ap- 
plications) can be found in [DBB06, PH10]. 


Radiance is influenced by the characteristics of materi- 
als with which electromagnetic radiation interacts by scat- 
tering, reflection and refraction. Reflection is the process by 
which electromagnetic flux (power), incident on a station- 
ary surface or medium, leaves that surface or medium from 
the incident side without change in frequency [NRH*92]. 
Our subject (the BRDF) is useful as a model for reflection. 
In general, and unless otherwise stated, we do not take into 
account subsurface scattering, refraction in the material, or 
participating media. We also assume all photons have a sin- 
gle wavelength A (fluorescence is not modelled), that reflec- 
tion happens in a short interval around an instant of time 
t (we omit phosphorescence), and that light is unpolarized. 
However some BRDFs models do take into account the ef- 
fect of some of these phenomena on the BRDF expression, 
and this review gives further details in those particular cases. 


In general, parameters t, x and À are implicit in our nota- 
tion, however in this section x is included explicitly for clar- 
ity. Regarding A, we exclude it from the notation for most 
BRDF models, except for the cases when the BRDF func- 
tion depends on it. 


When radiant power hits a smooth and locally planar in- 
terface between two media with different refraction indexes 
(the index of refraction of a medium is a value inversely pro- 
portional to light speed in that medium), a part of the incom- 
ing flux is reflected back to the incident medium and part is 
transmitted through the other medium (this part undergoes 
refraction, that is, a change in the orientation of the wave- 
front, governed by Law of refraction or Snell’s law). Well 
known Fresnel equations characterize the fraction of power 
which is reflected and transmitted, as a function of power, di- 
rection and polarization state of incident light [Hec02] (the 
term polarization state is related to the evolution in time of 
electric and magnetic fields vibration directions). 


We will use a simplified particle model for radiation, a 
model in which photons leaving (or reaching) x can be as- 
signed a direction u (a unit-length direction vector, with 
u € Qy, that is, it holds u- ny > 0, where nx is the unit-length 
vector perpendicular to the surface tangent plane at x). The 
reflection process causes reflected photons leaving x to be 
distributed in arbitrary outgoing directions. The distribution 
of directions for these photons depends on two factors: (a) 
the characteristics of the material at x and (b) the distribu- 
tion of directions for incoming photons at x. The BRDF de- 
pends only on the characteristics of the material. Note that 
Qx has an associated solid angle measure, which we will 
denote as 6. We will also use Gp to denote projected solid 
angle measure, that is, Op is the measure on Qx satisfying 
dop(u) := (u- ny) do(u), for any u € Qx. 


We can measure the total radiant power reaching a small 
area around point x per unit of time and per unit of area. 
This quantity is called the irradiance at x, written as E(x), 


its units are W-m~*. We may restrict this measure to the 
total summed power of all photons coming only from direc- 
tions in any given hemisphere region Rj C Qy. We will write 
this as ®;(x < R;). Symbol ®;, thus, stands for a measure in 
Qx. Total incident power ®;(x < R;) may cause a part of it to 
produce reflected photons (from x) in arbitrary directions. It 
is possible to consider which part of it is reflected to outgo- 
ing directions in another region Ro C Qx. We will write this 
power as F;-(x, Ro + Ri). Functions ®; and F have the same 
units as E£, and both can be considered as measure functions 
(in Qx and Qx x Qx, respectively). A detailed introduction to 
radiometric functions by using measure theory can be found 
in [Arv95a]. 


Both ®; and F, are absolutely continuous with respect to 
Op measure. This is because for any Rj, Ro C Qx, if Op (Ri) = 
O then necessarily ®;(x «+ R;) = 0, and if either op(Rj) = 
0 or Op(Ro) = 0 then F,(x,Ro + Ri) = 0. This continuity 
allows to use the Radon-Nikodym derivative (twice, one for 
each argument of F,) to define two new functions F! and F” 
as follows: 


dF, 

F} (x, Ro < Ri) = Jg (Ro + Ri) (1) 
p 
dF; 

Fy! (x,Wo < Ri) = Jo or We + Ri) (2) 
P 


Function F’ can be viewed as a measure function on its se- 
cond argument R; (it is the reflected power density in di- 
rection Wọ due to power incident from region Rj). Note 
that F” is the differential amount of power reflected at x 
to direction Wọ due to reflection caused by a differential 
amount of power incident at x from direction w;. Its units 
are Watt -meter~° - steradian’. 


Two additional functions L; and Lr can be defined by us- 
ing ®; and F/: 


d®; 
Li(x wi) := ie (x + wi) (3) 
L(x Wo) := Fi(x,Wo < Qx) (4) 


function L; is the density of power arriving at x from wj, 
and is called incident radiance. Function Ly is the density 
of power reflected from x to direction Wo, and is called re- 
flected radiance. Both quantities have units for W - m?. 
steradian—'. 


The Bidirectional Reflectance Distribution Function 
(BRDF), denoted by f+, is a function of two vectors W;, Wo € 
Qx, and can be defined formally as: 


Fy! (x, Wo = wi) 
Li(x 4+ w;) 


fr(X,Wo < Wi) := (5) 
this implies that we can view f- as the amount of reflected 
radiance at Wo due to a unit of irradiance at w;. Its units are 
steradian™ |. As a ratio of two non-negative power densities, 
fr is also non-negative and unbounded (although its integral 
is, see below). An important property is that the BRDF is 
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independent of L; (the distribution of incident radiance), and 
only depends on the characteristics of the material around 
surface point x. 


By using (2) and (5) it is possible to expand F! in (4) 
and express reflected radiance as a weighted integral sum 
of incident radiance, where the BRDF act as the weighting 
function: 


L(x Wo) := F(x, wo + Qx) (6) 


= l dF! (x,Wo < Wi) 
Q 


z 1 Fl! (x, Wo < wi)dop(wi) 


x 


= f fr(x, Wo < wi) Li(x < wi) dop(wi) 


The albedo (also called Bi-hemispherical reflectance) at 
point x is the ratio of total outgoing power to total incident 
power. It is written as p(x) (or simply p) and can be defined 
as: 


F(x, Qx — Qx) 

a 7 

PD oe) (7) 

Sometimes it is useful to consider the amount of total re- 

flected power in all directions due a unit of power coming 

from a single incident direction’ w;. This is the directional- 

hemispherical reflectance, pgp, and can be defined as fol- 
lows: 


pasew) = f fwo wi)dop(wo) ©) 


Note that both albedo and directional-hemispherical re- 
flectance depend only on the BRDF at x. 


In this paper, we consider each particular BRDF function 
as a formal model for a particular distribution of reflected ra- 
diance which happens in nature. This kind of models is usu- 
ally called physically plausible BRDFs, in contrast with non- 
physical BRDF models, which cannot exist in nature but can 
be used for applications where it is not necessary or desir- 
able to reproduce natural reflection. As stated, a physically 
plausible BRDF must be non-negative, but it also must obey 
these two additional properties: 


Symmetry: The Helmholtz Reciprocity Rule [Hel25] states 
that the BRDF is symmetric, that is: 


VuveQr: frx, uv) = f(x,y +u) 


this symmetry property leads us to use a notation where 
a double arrow is written between both vectors, thus we 
will write f-(x,Wo +> w;) instead of f(x, Wo < w;). It is 
important to state that this property does not hold always 
but just for certain polarization states of light [Hel25], 
[CP85]. However, in computer graphics and computer vi- 
sion literature, BRDF models are usually assumed to obey 
this property, although some researchers have explored 


t Li is an ideal Dirac-delta, zero everywhere except for w;. 
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the consequences of non-symmetric BRDFs in the imple- 
mentation of rendering systems [Vea97]. 

Energy Conservation: Energy conservation implies that 
no more radiant energy can be reflected from a point than 
the energy incident to that point, that is p(x) < 1. More- 
over, as this is true for any L;, it is in particular true for a 
Li in which all the energy comes from a single direction 
w;, and this implies that pgp (x + w;) < 1. From this and 
(8) we get the following inequality: 


VveEQ,: f fr(x,u v)dop(u) < 1 
uE, 


due to symmetry, this is also true if we integrate on the 
second variable instead. Many BRDF models include in 
their formulation a normalization factor —different for 
each one— necessary to correctly bound the resulting 
albedo. 


A BRDF is plausible if it is non-negative, energy conser- 
vative and reciprocal. BRDF functions can be categorized in 
two classes according to the presence (or absence) of rota- 
tional symmetry. Let us use the symbol Rota to denote a ro- 
tation transformation of œ radians (with ny as axis) for vec- 
tors in Qy. A material is called isotropic if f,(x,RotaWo © 
Rotqw;) is independent of a, and anisotropic when that does 
not hold. 


Total radiance L(x + wo) leaving x in direction wo is 
defined as the sum of emitted radiance Le(x —> Wo) (radi- 
ation produced in the material not due to reflection) and re- 
flected radiance L;(x > Wo). Moreover, incoming radiance 
at x from w; is equal to the outgoing radiance from another 
point y on the reverse direction (that is, L;(x + w;) = L(y > 
—w;)), because in empty space radiance is conserved along 
straight lines. Thus, by using (6), the value of the radiance 
function L at vacum can be written as: 


L(x > Wo) = Le(x > Wo) (9) 
J | frlX,Wo © w;) L(y > —w;) dop(w;) 
Qy 


This equation expresses radiance at one point as a weighted 
integral sum of radiance from other points in the scene and is 
called the rendering equation. It is an example of a variant of 
a Fredholm-type integral equation of the second kind, with 
fr as the kernel and L being the unknown function. It can 
be proved that a solution exists and it is unique, provided f- 
obeys conservation of energy (see [ATS94], section 5.3, note 
that pan must be strictly less than one). 


As has been stated, the main goal in Global Illumination 
is computation of values of the radiance function L. How- 
ever, this is a complex calculation because in Eq. (9) func- 
tion L appears on both sides, that is, it can be interpreted as 
a recursive definition of L. Radiance is usually computed in 
Global Illumination by using Finite Elements methods and 
Monte Carlo (MC) integration techniques, because it is of- 
ten impossible to obtain analytic expressions for L [DBB06]. 
When using Monte-Carlo algorithms for rendering, an im- 
portant step is BRDF sampling, which is described in more 
detail in section 5 of this paper. 
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2.1. BRDF parametrization 


In this section we introduce several alternative BRDF param- 
eters domains or parametrizations. Original BRDF defini- 
tion in (5) clearly implies that the domain is Qy x Qx. How- 
ever each BRDF model can be expressed by using different 
formulae depending on how vectors in Qx are represented. 
This leads to different formulae complexity, and different 
accuracy when fitting data based on measurements. More- 
over, it yields also different time complexity and numerical 
accuracy in BRDF-related computations. 


Parametrizations make use of various particular direction 
vectors in Qy, which we introduce here: 


e Tangent vector (t), this vector is perpendicular to ny and 
thus tangent to the surface at x. It is used to build a lo- 
cal reference system. BRDF expressions can be refered 
to. Isotropic BRDFs may use any tangent vector, however 
for anisotropic BRDFs this vector must be considered as 
an external parameter which fixes the orientation of the 
BRDF (with respect to rotations around nx). 

e Bitangent vector (b), defined as b := nx x t. 

e Halfway vector (h), defined as h := (w; + Wo) /||wi + woll 

e Vector b’, defined as ny x h/||nx x h||. In the case ny = h, 
vector b’ is defined equal to b. 

e Vector t’, defined as b’ x h 


Parametrizations are usually referred to local reference 
frames. By using the vectors introduced above, we can de- 
fine reference frame Re fn, whose X,Y and Z axes are aligned 
with t, b and nx, respectively. We can also define reference 
frame Re fp, whose Z axis is aligned with the halfway vector 
h, Y axis is aligned with b’, and X axis is aligned with t’. 
These reference frames allow us to define various possible 
parameters domains in Qx. 


Global cartesian coordinates: in these cases, vectors w; 
and wo are expressed as 3-tuples containing world coordi- 
nates, that is, coordinates relative to world reference sys- 
tem used to express objects geometry in the scene model. 

Local cartesian coordinates: now, w; = (xi,Yi,zi) and 
Wo = (Xo, Yo; Zo) are again expressed as 3-tuple, however 
these tuples contain coordinates relative to local reference 
frame Re fn. Thus, the following equalities hold: 


Xj i= Wi t yii=wi-b Zi i= Wi: Ny 


10 
Xot=Wo-'t Yo:=Wo:b Zo:= Woy (10) 


Spherical coordinates: in this case, vectors are expressed 
as 2-tuples: w; = (Ọ;,9;) and Wo = (@o, 80). Each pair 
contains the azimuthal angle and normal angle, respec- 
tively (see Figure 2.1, left). These angles can be defined 
from local coordinates (as defined in (10)), as follows: 


Ọ; := atan2(y;/x;) ©; := arccos(z;) 
Qo := atan2(yo/x0) 90 := arccos(zo) 


(1) 


Spherical coordinates for isotropic BRDFs: these 
BRDFs depend only on the difference between Q; and Qo, 
but not on the two individual values. Thus it is possible to 
express the BRDF just in terms of three values, namely 
0;,0, and Q4, where Oa := Qo — Q; 


Half angle- difference parametrization: this parametriza- 
tion was introduced by Rusinkiewicz [Rus98], who 
showed it is well suited for functional approximation of 
reflectance data based on measurements (see section 4.2). 
The BRDF is expressed in terms of the spherical coordi- 
nates of h relative to Refn, which are written as (@p,9,), 
and the spherical coordinates of w; relative to Re fp, writ- 
ten as (@4,9,) respectively (see Figure 2.1, right). Thus, 
these angles can be defined as: 


Qh = atan2(yp/xp) On := arccos(zp) (12) 
Og := atan2(yq/Xxqa) 07 := arccos(z,) 

where: 
xX,i=h-t yn i=h-b Zh := h -ny (13) 
xa:= Wit ya:=Wi b z:=wjeh 


Note that isotropic BRDFs are independent of Ẹ}, thus 
these BRDFs can be expressed just in terms of 0,04 and 


a. 


Figure 1: Spherical coordinates, including azimuthal angle 
(®), and elevation angle (8), for vectors in Qx. Coordinates 
can be relative to the reference system Refn, with Z axis 
aligned with nx (left) or relative to reference system Re fh, 
with Z aligned with h (right). Red arc includes vectors with 
o=0. 


2.2. General classification 


In general BRDF models can be classified into one or two 
of these categories: empirical, theoretical and experimental. 
Moreover, a BRDF can be based on a previous one, extend- 
ing their capabilities and the type of material that it is able to 
model. Figure 2 represents a graphical classification of the 
BRDFs which are introduced in section 3. 


Empirical 
Their main aim is to provide a simple formulation specif- 
ically designed to mimic a kind of reflection. Conse- 
quently, we get a fast computational model adjustable by 
parameters, but without considering the physics behind it. 

Theoretical 
These models try to accurately simulate light scattering 
by using physics laws. They usually lead to complex ex- 
pression and high computational effort, thus they are not 
normally employed in rendering systems. 

Experimental 
The BRDF can be acquired using a gonioreflectome- 
ter [War92, GA97, GA99, HLW06] which mechanically 
varies light source and sensor positions. This process 
could take hours and usually data is limited by some 
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angular resolution. Other techniques use digital cameras 
to acquire many BRDF samples with a single photo- 
graph [DGNK99, MPBM03]. No much densely acquired 
data is readily available. 


We have the power of using the reflectance data directly 
in the rendering process as [MPBMO03] obtaining great real- 
ism in the visual results. The other option is to approximate 
the reflectance data with an analytical expression. We get a 
formula which is simple, fast, accurate and adjustable by pa- 
rameters. Some reflectance models (empirical or theoretical) 
are being designed for this use in mind. In these cases the 
parameters are not necessarily intuitive as they are meant to 
be set algorithmically to fit measurement data. Additionally 
the fitting of data presents several problems [KC08]: (1) the 
measurement error introduced in data could make it difficult 
to find a better fit for your model, (2) the solution to the fit- 
ting process (global minimum) is not guaranteed, (3) is not 
easy to know what is the best objective function, and what 
are the initial values to use, and (4) in some cases, you need 
to find a suitable number of lobes (it is recommended three 
or four in the Lafortune BRDF [NDM05)), at the expense of 
the non-linear minimizer becomes unstable and fails to find 
the minimum objective function. 


3. A review of BRDF models 


In this section we make a review of some of the models 
proposed in the literature, following the classification sug- 
gested in Figure 2. We introduce the desirable properties for 
BRDF models, and then we analyze if a particular model 
fulfill them and in what degree. Table 1 summarizes some of 
characteristics of the reflectance models reviewed in this pa- 
per, including their relative computational cost considering 
that we have implemented 16 of them. 


It is understood in Computer Graphics that reflectance 
models should exhibit a set of desirable properties [Shi96] 
to make them realistic and reliable at the same time: 


1. Physically plausible: a function that obeys non- 
negativity, reciprocity and the law of energy conserva- 
tion. A BRDF with this property can be used safely in 
a rendering system, avoiding situations where energy is 
created wrong. 

2. Expressive: the model is adjustable by parameters. These 
should be intuitive to set. 

3. Usable: capable of representing as many different mate- 
rials as possible, including those acquired with a device. 

4. Realistic: close to actual BRDF functions which can be 
found in nature, showing directional diffuse and Fresnel 
behaviors. 

5. Efficient: must lead to efficient implementations on 
Global Illumination rendering algorithms and thus quick 
to evaluate and appropriate for importance sampling of 
Monte-Carlo integration. 

6. Accurate: Reflection components (diffuse, directional- 
diffuse, specular) should be represented by the model it- 
self, avoiding ideal simplifications. 
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A 
2 3 g 
o $ 8 2 
z 2 m 6 oo z = 
2 ¢ 6 2 & g Š 
Models £ A È <x & 4 = 
Ideal Specular * * v v * x perfect specular 
Ideal Diffuse v v x perfect diffuse 
Minnaert M a v v v 5.35x Moon surf. 
Torrance-Sparrow * M * * M E Tough surf. 
Beard-Maxwell * v * M v 397x painted surf. 
Blinn-Phong M v v M * 9.18x rough surf. 
Cook- Torrance * * * M v 16.9x metal, plastic 
Kajiya * M * * M a metal, plastic 
Poulin-Fournier * M v * v 67x clothes’ 
Strauss v ve * v v T4.B8x metal, plastic 
He et al. * * * v v 120x metal 
Ward v v v OX woo 
Westin * vee * * v vee metal 
Lewis v * v v * T0.73x mats 
Schlick Y * * * v 26.95x heterogeneous 
Hanrahan * vee * v v vee human skin 
Oren-Nayar * * v v * T0.98x matte, dirty. 
Neumann y * y * * vee metal, plastic 
Lafortune v * vy * * 5.43x rough surf. 
Coupled * * * v * 17.65x polished surf. 
Ashikhmin-Shirley * v * * * 79.6x polished surf. 
Granier-Heidrich * aia * v v tee old-dirty metal 


Table 1: Brief summary of the properties exhibited by the reviewed 
BRDFs. Legend: (X) if the BRDF has this property; (¥) if the BRDF 
does not; (:-- ) unknown value. 


3.1. Physical-based Reflectance Models 
Ideal Reflection 


In case of ideal specular reflection, light incoming from a 
given direction is reflected in a single direction following 
the law of reflection. The BRDF in this case is a delta dirac 
distribution, 6, giving always zero, except for the reflection 
direction r. This reduces the radiance computation consider- 
ably since L;(x > Wo) = Li(x + rw,)- 


fr(Wo, Wi) = Ps(Wi) O(Wo, wi) 


where ps is the specular reflectance at the point and 6 is de- 
fined as follows: 


def co ifv=r 
d(ru, v) a { 0 7 


For this BRDF, its hemispherical integral is 1 when the inci- 
dent and the outgoing direction are at the same angle relative 
to the surface normal and 0 otherwise. 


A diffuse surface has a BRDF that has the same value 
for all incident and outgoing directions. This substantially 
reduces the computations and thus it is commonly used to 
model diffuse surfaces as it is physically plausible, even 
though there are no pure diffuse materials in the real world. 
This BRDF is expressed as: 


Aj = Pa 
fr(Wo, Wi) — T 


where pq is the diffuse reflection? . It includes the value 7 for 


+ The fraction of energy reflected w.r.t. the total incident energy, 


that is Z-O 22o) 


PEER) by using the terminology introduced in section 2. 
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Figure 2: A graphical classification of the BRDFs cited in this paper. Some BRDFs are built on previous ones. 


normalization, so is the result of integrate cos(v) in the hemi- 
sphere of directions to compute the reflectivity or albedo. 


The Torrance-Sparrow BRDF 


This BRDF is one of the most complete physical reflection 
models for isotropic materials. It is considered precursor to 
other models and its formulation has been validated by a 
ray-tracing-like simulation [MWM* 98]. It considers polar- 
ized light and is used for rough surfaces [TS66, TS67]. The 
roughness is modelled using microscopic concavities in V- 
form of equal length called microfacets. Their orientation is 
random and their distribution is controlled by parameters, so 
it is possible to simulate different degrees of roughness. The 
complete BRDF function is as follows: 


ka ks 


= aw G(wo, wi) 


fr(Wo, Wi) = 


We consider the previous terms separately: 


e The microfacets distribution D : Q — R with fo D(h) = 1 
gives a distribution of the normals of the microfacets that 
are aligned relative to vector h and it is parameterized by 


m. Many authors use the Gaussian distribution function 
but the Beckmann distribution [BS63] is also common: 


a cos(8)* — 1 
m? cos(6)* E m? cos(8)? 


D(h) = 


The Fresnel factor F (wọ) € [0,1] gives the fraction of 
light that is reflected from the entire surface. Its computa- 
tion is a linear combination of the coefficient for perpen- 
dicular light polarization and the coefficient for parallel 
light polarization. It is quite usual to use the Schlick ap- 
proximation of this computation [Sch94b]. 

Geometric attenuation factor G(wo, w;) € [0, 1] expresses, 
for two directions, the ratio of light that is not occluded 
by the surface due to masking or shadowing. This model 
takes these effects into account by evaluating the follow- 
ing formula: 


G(wo,Wi) = min {1,2M Wo a 


(wo:h) °  (wo:h) 
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This microfacet model was the basis for many other works 
who offered variations on the calculation of the functions D, 
F and G. Figure 3 shows the vector system involved in the lo- 
cal reflection at a given surface point, and will be referenced 
in subsequent reflectance models. 


Figure 3: Angles relative to incident and reflected vectors in a local 
coordinate system. 


The Beard-Maxwell BRDF 


With the description of this BRDF [MBW*73] one can 
simulate, based on physical properties, a very concrete type 
of material: painted surfaces. The final model considers the 
combined contributions of a local specular reflection on the 
first layer (whose normal is n), named f; sup, and a volume- 
tric reflection approximation which comes from the reflec- 
tion of light in the interior layers, simulating subsurface 
scattering and noted by f, vol- 


Sr (x, A, Wo, Wi) = Frsup (A, Wo, Wi) + frvol (A, Wo, Wi) 


The superficial component computes a specular reflection 
that takes place on the first layer of the surface. The first 
layer is represented by microfacets whose normals are ori- 
ented following a statistic distribution D, as the used by 
Torrance-Sparrow. It is also governed by the Fresnel term 
F for dielectrics (using K = 0 and n = 1.65): 


F(B) _fr(h) cos*(h) 
F(0) cos(wo) cos(w;) 


frsup (Wo, Wi) = SO(wo,h,T, v0) 
F (0) D(h) 
4 cos(Wo) cos(w;) 


frh) = 


As with Torrance-Sparrow [TS66, TS67], the surface is 
modelled with microfacets so the same shadowing and 
masking (here is called obscuration) could occur. Though 
the light is blocked in the same form, the formulation that 
models this effect is not the same. The term SO summarizes 
both situations, and is dependent of two parameters T and v 
that take their values from measured data: 


ee e7 1 
SO(wo,h,t,v) = s n 
bp E oe 


v 


The volumetric component is a simple model of the sub- 
surface scattering of the light, where light is supposed to be 
depolarized: 


2 pv f(B) s(n) 


a w. w) = 
Frvor rl 0; i) cos(wo) cos(w;) 


where py is the measured reflectance of the surface when 
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Ow, = Ow; = 0 and f(B) = g(h) = 1 is used. For a com- 
plete description of the f and g functions please refer to the 
original paper [MBW*73]. 


Westlund and Meyer published a database of 400 mea- 
sured materials which fitted with a modified Beard-Maxwell 
BRDF named NEF-BM [WM02]. This modification simpli- 
fies the original formulation by dropping the parameter op 
and the f and g functions. 


The Cook-Torrance BRDF 


Based on geometrical optics theory and in previous stu- 
dies [BS63, TS67], Cook & Torrance developed a model 
that introduced new ideas: only those microfacets oriented 
toward h vector contribute to the final reflection. This 
model [CT82] also introduces a new type of material in the 
field of rendering, differentiating between metallic and non- 
metallic surfaces. 


Reflection is described using a combination of the diffuse 
and specular parts with parameters ky and ks controlling the 
fraction of energy that is diffusely or specularly reflected. 
This parameters are considered a characteristic of the mate- 
rial. The diffuse term is represented by a classical Lamber- 
tian reflection. The specular term is the compound of the al- 
ready known functions F, D and G of the Torrance-Sparrow 
model [TS67]. In this case, the normalized halfway vector 
simulates an imaginary surface which behaves as a mirror, 
and the Fresnel term represents the reflection of polished sur- 
faces. 


One of the major contributions of this work [CT82] was 
the development of an optimized formulation of Fresnel, 
making this model more applicable to Computer Graphics. 
In the original formulation, the energy reflected depends on 
the index of refraction n, and the absorption coefficient «K. 
For some materials both values are unknown, but could be 
approximated by assuming that the value of the reflectance 
at normal direction is constant. This facilitates to solve the 
coefficients dependent on the wavelength: a) n is calculated 
assuming that F(}) = 1 and « = 0; b) x is calculated by 
knowing F (0) and having n = 1. Finally the simplified Fres- 
nel expression is: 


F(e) = 1 (b= cy" fis 


2 (b+c) [e(b—c) +1] 


[c(b +e) — 1]? ) l 


where b? = n? +c? —1andc = cos(8). 


Some surfaces have two or more scales of roughness, con- 
trolled by the slope m (with m € IR), and can be modelled 
by using two or more distribution functions [BS63]. In such 
cases, D is expressed as a weighted sum of the distribution 
functions, each with a different value of m parameter. When 
m is small, the microfacet normals distribute mainly towards 
the reflection direction. As m value grows, the distribution 
expands and becomes more uniform. Typically, D is a Gaus- 
sian distribution: 
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The geometrical attenuation factor G accounts for the 
shadowing and masking of one microfacet by the other side 
of the microfacet. Its definition is identical to the Torrance- 
Sparrow model (see section 3.1, page 6). Finally, the specu- 
lar component is: 


F(B) D(h) G(wo, wi) 
T cos(Wo) cos(w;) 


frs(Wo, Wi) = 


The main disadvantage of this model, compared to other 
BRDF representations, is that the parameters are not intui- 
tively set. The user needs to experiment with the values until 
a good result is found. Another drawback of the model is 
that the function is not plausible because for some angles 
this BRDF does not obey the energy conservation law. 


The Kajiya BRDF 


The Kajiya model [Kaj85] implements an anisotropic 
method that estimates the analytical form of the reflected in- 
tensity light. It is based on Kirchoff’s approximation and a 
stationary base method for the approximation of the radi- 
ance equation (Eq. 9), and is considered a model with high 
computational costs as it is noted in table 1. 


The model considers a simplified rough surface that uses 
the nearest tangent plane. In this work, Kajiya calculates and 
stores the reflectance in a table each time a surface illumi- 
nation is evaluated. Later, values from the table are used to 
perform linear interpolation. Though he uses a novel numeri- 
cal technique that explores unpolarized light’s properties and 
the Fresnel term, this model shows a lack of many important 
inter-reflection phenomena. Moreover, conservation of the 
energy is not guaranteed. 


The Poulin-Fournier BRDF 


Poulin and Fournier [PF90] introduced a model for 
anisotropic materials that considered reflection and refrac- 
tion of the incident light. Anisotropy is represented by a set 
of parallel cylinders packed tightly together such as those 
used in [MH84], whose optical properties —and thus the 
degree of anisotropy— are described by varying the height 
and distance of these cylinders. The surface is modelled by 
a distribution of cylinders that are added or subtracted from 
the base surface. The phenomenon of reflection or refraction 
happens at the longitudinal cut of the cylinder, more specifi- 
cally, at the binormal plane BN where b = nx t, so we 
assume the local system is formed by n, t and b vectors (see 
Figure 3 right in page 7). 


To control the degree of anisotropy, this model uses 
two parameters: distance between cylinders d € [0,co], and 
height from the base surface h € [0, 1]. Both are unit less and 
appear in Figure 4. When d = 0, this model is the Torrance- 
Sparrow BRDF. When d = 2, the normal of the cylinders 
differentiate in 1 and we get the maximum anisotropy. How- 
ever, if d > 2 an imaginary base surface with a constant nor- 
mal vector appears, and takes a role in the reflection of the 
light. The base surface is modeled with a classic combina- 
tion of diffuse reflected intensity cos(w;) and the cos(h)” 


term of the Blinn model to consider the light reflected specu- 
larly. The height of the base surface is then controlled with 
the h parameter. Increasing h causes the decrease of normal 
variation, thus softening the anisotropy effect. 


Base 
Surface 
h 

\ 

y 


bg eenma i <--- -e 
d radius=1 


Figure 4: Anisotropy is controlled with the distance and height pa- 
rameters. 


The amount of reflected light is calculated by determin- 
ing the visible and illuminated portion the cylinder on the 
incident light, as well as reflected by the visible part of the 
base surface when it intervenes. To make these calculations 
is necessary to know the length of the projection of the cylin- 
der on the base, visible and illuminated part of the cylinder, 
as well as the visible and illuminated part of the base sur- 
face. The model also considers the shadowing and masking 
effects of light. You need to know the angles where the cylin- 
der begins to hide himself from the direction of observation 
Osh, or by other cylinders. Similarly angles are calculated at 
which cylinder starts to block light toward the viewer, due to 
himself Oss, or by those cylinders at its side. 


The calculation for the refraction event is the same with a 
simple exception: the cylinders are negative or modelled in 
reverse. In Figure 6 third column two examples of the plot of 
this function are shown. This model is not energy conserving 
and its visual results depends upon knowing how to set the 
parameters correctly, as for the rest of the BRDF models. 


The He-Torrance-Sillion-Greenberg BRDF 


One of the most complete and complex physical BRDF mod- 
els is the one of He et al. [HTSG91]. It considers impor- 
tant phenomena associated with the wave-like nature of light 
such as diffraction and interference. Many parameters con- 
tribute to the representation of isotropic materials. Some are 
geometrical, such as angles, and are illustrated in Figure 3, 
and others are physical, for example the wavelength A of the 
incident beam, the index of refraction y (used in the Fresnel 
term calculation) and parameters © and T to characterize the 
roughness. 


This BRDF starts by distinguishing three components in 
the reflection event: the uniform diffuse reflection, the per- 
fect specular, and the directional diffuse reflection com- 
ponent. The authors describe an implementation based on 
spherical harmonics decomposition (see section 4.4) and the 
Kirchoff’s theory for the rough surface analysis. They are 
used both to approximate the reflectance distribution func- 
tions pa and fe and directional intensity distributions for 
illuminated surfaces. The perfect specular component f,” is 
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represented by the Fresnel equations. By summing the ex- 
pression of the three components we get the complete BRDF 
equation: 


[Fw JP? exp(—8(6.A)) S(1) 


fr(Wo,wi) =a(A)+ cos(wi) oF 
2 a m exp(—o(o,d 242 
area ` 16% “Lin=1 s atte a AR (“in 


where a is the ambient wavelength function, g is the effective 
surface roughness, © is a parameter for the amplitude of the 
roughness, S$ is a shadow function, T is a parameter defining 
the autocorrelation length, 6 is a delta dirac function, F is a 
function derived from the Fresnel function, b is the bisected 
unit vector, p is the incident polarization state vector and w 
is the wave vector change. 


Detailed mathematical construction is included 
in [HTSG91], experimental analysis can be found 
in [NDM05]. Although at present is not an issue to 
consider, at the time of this model the computational cost 
was high and required optimizations based on precalculus 
and search tables [HHP*92]. In attempts to make it more 
usable, simplifications have been made in the formula- 
tion [WAT92] and spherical wavelets has been used to 
represent this BRDF. This idea is detailed in section 4.5.1. 
The He et al. BRDF is theoretical, not reciprocal and since it 
has no associated direct sampling method (as table 1 shows) 
is not usually integrated into Monte-Carlo based systems. 


The Oren-Nayar BRDF 


The Oren & Nayar model [ON94,ON95] is an improvement 
on the classical Lambertian interpretation for a diffuse ma- 
terial. This model is able to explain the view dependence ap- 
pearance of the matte surfaces with geometric optics, and for 
this reason the BRDF is classified as physical or theoretical. 


In this case the microfacet or groove distribution is purely 
diffuse. In addition to masking and shadowing light block- 
ing events, the reflections between microfacets are also con- 
sidered, but they have a limited number of bounces within 
each pair of microfacets. The orientation of the grooves is 
distributed according to a Gaussian function with a standard 
deviation for the angles Œm (in radians), which is a mea- 
sure of the roughness of the surface. It ranges from 0 to 1. 
The dirty materials tend to more reflect the light to backward 
of the light(retroreflection). This characteristic is controlled 
with parameter On. The larger Om, the more retro-reflective 
is the material. 


Given an incident direction, the total reflected radiance is 
the integral of grooves’ reflectance values —the projected 
reflected radiance— and the radiance due to the bounces be- 
tween them. Oren & Nayar found hard to solve the result- 
ing equations analytically, so they presented a functional ap- 
proximation: 

p 


fr(Wo, Wi) = p (A+B max(0,cos(w; — Ow, )) sin(a) tan(b)) 


where a = max(@w, , Ow, ), b = min(Ow, , Ow, ), and: 
o, 7 o, 
0.5 F033 B= 0.45 27009 


m 


A=1 
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In the case of Om = 0, roughness is null and the previous 
expression is equivalent to the Lambert Law (also A = 1, 
B= 0). Application of trigonometry transformations can 
substantially improve the implementation of this BRDF. 


The Coupled BRDF 


Shirley et al. [SHSL97] introduced a model for the represen- 
tation of matte surfaces. It is called a coupled model because 
it properly combines the diffuse and specular reflections. It 
is plausible, physically based and suited for Monte-Carlo 
rendering algorithms. The model for specular reflection as- 
sumes a polished planar surface of dielectrics in touch with 
the air. The matte component is considered to be almost con- 
stant and isotropic, so both are combined as follows: 


Sr(Wo,Wi) = F(8w,) +kRm f(Ow,) f(8w;) 


where k is a constant for normalization and Rm is the re- 
flectance of the matte component given as a parameter of the 
model. Finally f(O) is approximated with the Schlick ap- 
proximation of the Fresnel term [Sch94a] as follows: 


f(8) x (1 — (1—cos(8))”) 


Restrictions are imposed to guarantee energy conservation 
as if the Fresnel term at grazing angles tends to one, the dif- 
fuse component should tend to zero: 


NIA 


Ro + 27k f (Ow; )do(w;) =1 
0 


allowing the authors to solve k: 
= 21 
~ 20n(1— Ro) 


where Rọ is another parameter that varies between 0.03 and 
0.06. We arrive at: 


fr(Wo,wi) = [Ro + (1 —cos(wo))° (1 — Ro)] frs(Wo, Wi) + 
kRm [1 — (1 —cos(wo))°][1 — (1 — cos(w;))5] 


The Ashikhmin-Shirley BRDF 


Ashikhmin & Shirley [ASOO] developed an anisotropic mi- 
crofacet model that uses simple and intuitive parameters. Ba- 
sically you can generate any BRDF from a distribution of 
microfacets which is generally expressive. They also pro- 
vided a method of sampling based on the BRDF function. 
Two parameters define the anisotropy by controlling the 
halfway vector distribution of the microfacets. These para- 
meters are related to the axes of an ellipse: ex and ey that 
orientate h along the X and Y axes respectively. An expo- 
nent value e is calculated with respect to the ellipse radius 
and the @ angle of orientation. The distribution function is as 
follows: 


D(h) = y/(ex + 1) (ey + 1) (h-m) 0 On) +60" (n) 


Applying trigonometry rules, a simplification of the expo- 
nent of D(h) is achieved. Note that the following expression 
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are undetermined when zp = 1. 


sin” (On) = iz 


2 


cos?(Mn) = phy 


In [APSOO] the authors complement the previous BRDF 
model with one for diffuse surfaces covered with a polished 
or painted layer. For example a wooden table with lacquer 
paint. At normal incidence it is diffuse but at the horizon- 
tal view it is almost specular. The Fresnel equations describe 
this very phenomenon, and it is modelled with a superposi- 
tion of two layers: one specular and one diffuse. 


The BRDF is expressed as a weighted sum of the specu- 
lar term and the diffuse one: f = frs + fra. For the specu- 
lar reflection, the distribution over h is used as well as the 
Schlick’s approximation of the Fresnel term [Sch94b]: 


D(h) F(wi) 
87 (h- wo) max(cos(Wo),cos(w;)) 


frs(Wo, Wi) = 


The diffuse term guarantees the reciprocity and energy 
conservation properties, so the BRDF is plausible. 
28kq 


Sid (Wo, Wi) = 33n (1 —ps) A(wo) A(wi) 


nia) (1 (1 =w’) 


Two examples of the function plot are shown in Figure 6 
first column. Though this BRDF is not symmetric is one of 
the most complete reflectance function as it also includes 
a sampling method for the generation of random reflected 
directions. 


where 


Granier-Heidrich BRDF 


The Granier & Heidrich BRDF [GH03] is a reflectance 
model specifically conceived for an easy to evaluate repre- 
sentation of old metals, or dirty surfaces covered with an oily 
transparent exterior surface. This two layers model uses the 
wavelength of the incident light to simulate two phenomena: 


1. A change of phase between the beam of light that is re- 
flected and the beam of the incident light that is transmit- 
ted to the other layer occurs. 

2. The index of refraction is wavelength dependent, n (À) 
and produce the dispersion of light in different colors 
(much like a prism). This is obtained with a model of non- 
parallel layers, with a different alignment of the normal 
of the outer surface and the normal of the inner surface. 


Though the model is à dependent, an initial approxima- 
tion is carried out using the RGB color model with values 
R = 645 nm, G = 525 nm, and B = 445 nm. Another sim- 
plification is associated with the outer layer using N; cor- 
responding to the vacuum, with nı < n2. This assumption 
implies that transmission always happens. The final result is 
a BRDF that combines the reflection term R and the trans- 
mission term T. R is a cosine weighted by the reflectivity 


—the Fresnel term using the Schlick approximation— due 
to direct illumination. For the transmission term, the reader 
is referred to the original paper [GH03], where a hardware 
acceleration technique that uses textures, is also suggested 
in order to obtain low cost results. 


3.2. Empirical Reflectance Models 
The Minnaert BRDF 


One of the first reflectance models was given by Min- 
naert [Min41] to model the lunar surface reflectance. This 
model can be applied not only to the moon, but to shade 
an object where we would see a darkening of color near the 
edges. The model is controlled using two parameters: pg and 
an exponent responsible for darkening k. This shape parame- 
ter usually varies between 0 and 2. When k = 1, this function 
is equivalent to the Lambertian function. The analytical ex- 
pression for this BRDF is: 


frlWo,Wi) = Pa (cos(Wo)cos(w;))*~! 


The Phong BRDF 


Phong [Pho75] is still a very popular BRDF model and it was 
the first description for non-Lambertian surfaces. Is a well- 
known class of BRDF models based on cosine-lobes where 
the term lobe here refers to the shape of a classic Phong-like 
BRDF as is shown in Figure 6 second column. Basically, it is 
an empirical model which obeys neither energy conservation 
nor reciprocity, but its simplicity has made it one of the most 
used in Computer Graphics. It depends on an & angle, which 
is calculated after reflecting the incident direction. Essen- 
tially, this model is a simplification of the Torrance-Sparrow 
one, where the G and F factors are dropped and the distribu- 
tion D is reduced to: 


D(Wo, Wi) = (Wo «Kw; )” 

where the parameter n € [0,00[ characterizes the shape of 
the specular highlight (from O or dull to more glossy sur- 
face). Furthermore, this model could be faster with an opti- 
mization of the exponential operator, like the one given by 
Schlick [Sch94a]: 


cos() 
n — ncos(a) + cos(a) 


cos" (a) ~% 


The second and third parameters are the specular ks and dif- 
fuse ky constants, both taking values in [0, 1]. They allow a 
linear combination of the specular lobe with the Lambertian- 
diffuse BRDF. Precisely the lack of physical significance and 
the fact of combining the two components linearly, can cre- 
ate situations in which energy is not retained. 


fr(Wo,Wi) = ks (Wo: rw)” 
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The Blinn BRDF 


The Blinn BRDF, also known as the Blinn-Phong reflection 
model, is the standard lighting model used in both the Di- 
rectX and OpenGL rendering pipelines. Blinn [Bli77] noted 
that normally the direction that gave the higher reflection 
correspond to the halfway vector h aligned with the normal. 
This vector is then used in conjunction with the normal to 
get the specular highlight. 


D(h) = (n-h)” = cos(8)" 


In this way, Phong BRDF is modified replacing (u- r) by 
(n-h), allowing less calculations, because it is not necessary 
to find the reflection vector. Two examples of the evaluation 
of this function is shown in Figure 6 second column. 


fr(Wo,Wi) = ks(n-h)” with kg tks = 1 


The Lewis BRDF 


Reviewing the previous models, Lewis [Lew94] made a 
study of their properties. Models like Phong or Cook- 
Torrance did neither exhibit energy conservation nor reci- 
procity, and so Lewis elaborated a plausible version of them. 
With reference to microfacets based surface modelling, the 
cause of plausibility is the normal distribution D, so it should 
be normalized in order that the BRDF can be used safely in 
any lighting algorithm. The total value of the projected mi- 
crofacet’s areas should be same as the differential area of 
the considered surface. If we assume energy conservation, 
the following should occur: 


[ D(h) cos(h) do(h) = 1 
Q 


In the case of the Phong distribution D(u, v) = (u- ry)”, the 
normalization factor isc = cae If we apply this normaliza- 
tion to the Blinn BRDF model, then energy conservation is 


guaranteed. We name Lewis BRDF to the modified version: 
n+2 ñ 
n-h 
zz (ah) 


We note here that a normalization to the Phong BRDF was 
given also by Lafortune et al. [LFTG97]. 


fr(Wo, Wi) = ks 


The Neumann-Neumann BRDF 


The Neumann-Neumann BRDF model [NN96] explores a 
plausible function easy to integrate in an Monte-Carlo based 
rendering algorithm. For the description of metal, plastics, 
ceramic, or retroreflective material, and even the anisotropic 
reflectance, the authors built these BRDFs over a base 
model. They presented two key ideas. The first, was that the 
base BRDF uses the unit disc domain D? and thus the Op 
measure, so a cosine term is implicit. In the second one, the 
projection P of the reflection vector is the origin of a disc 
C whose radius is controlled by a parameter r (r 4 0). Fig- 
ure 5 shows the geometry previously described. The BRDF 
function is a constant value defined as: 

fr(Wo, Wi) = { az POE ECP Naat) 


in other case. 
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Figure 5: The Neumann-Neumann BRDF is based on a disc of 
variable radius defined over D? domain. 


The surface may vary from perfect specular r = 0 + £, to 
perfect diffuse r = 2 and the model allows to simple calcu- 
lation of the directional-hemispherical reflectance: 


P(Wo) = r-Area(C(P(tw,),7)). 


Using the base model it is easy to build other BRDFs. For ex- 
ample a retroreflective one would simply change the h of the 
distance function by g = Wo — wj. As a consequence, light 
will reflect predominantly in the incident direction. Another 
example is a BRDF for anisotropic surfaces which change 
the circle C by an ellipse E whose axes are A and B, and the 
distance function by the normal elliptic curve. 


ak if P(wi) € E(P(tw,),4,b), b < B/A 
vers if P(wi) € E(P(tw,),4,5), B/A <b<2 


franiso (Wo, Wi) = 
0 in other case. 


The Strauss BRDF 


The Strauss BRDF [Str90] is an empirical model that cha- 
racterized metals, and other glossy surfaces that exhibit off- 
specular peaks. Strauss noted on previous BRDF models 
how difficult it was to select the correct values for those 
BRDFs when wishing to mimic the appearance of a mate- 
rial. His model is based on few and very intuitive parameters 
that represent different types of material properties that any 
user should be able to set without problems. 


e Chromatic color C or surface color (in [0, 1]). 

e Specular color C; is the calculated (not a parameter) shiny 
color. 

e Smoothness: when s = 0 we have the perfect diffuse sur- 
face, and when s = 1, the surface behaves like a mirror. 
This value changes the ratio from diffuse to specular and 
also the shape of the specular highlight. 

e Metalness: when m = | the surface is full metallic, and 
when m = 0 we have the absence of this property. It 
changes the diffuse aspect of the resulting image and also 
it is involved in the computation of the specular color. 

e Transparency: when t = 0 you get a totally opaque sur- 
face, when t = 1 all light not reflected will be transmitted 
through the media. The transmittance value used is the 
value of t at normal incidence direction. 

e Index of refraction n: employed to compute the direction 
of the transmitted beam of light. 
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The Strauss reflectance model uses a linear combination of 
diffuse and specular components controlled by two coeffi- 
cients ky and ks. 


fr(wo,Wi) = = ka fra(u,v) +ks frs(Wo, Wi) 
ka cos(w;)d ra C + ks rs Cs 


In order to obtain the values needed, the model performs 
some simple operations: 


ra = (1—s°)(1—1) d = (1—ms) 

rn = (l—t)-r4 e=74 

rj = min[l,rn+ (ra+k;)j| kj =0.1 

rs = [-(h-wi)]°rj j = F(a) G(a) G(8) 


The meaning of r; is similar to the Fresnel term together 
with the geometric term. It uses two constant values kp = 
1.12 and kg = 1.01, as well as a cosine term x € [0,1]. The 
model precomputes F and G, keeping its simplicity. It also 
assumes that for metals, the color of the reflection should not 
be the same as the color of the surface, and thus: 


if non— metal 


Gat 4 
aaa { Cj +m(1—F(a))(C—C,) if metal 


where C4 is the equivalent to white specular color that is 
acquired at 1/2 incidence angle. 


3.3. Experimental Reflectance Models 
The Ward BRDF 


Previous BRDF models that account for anisotropy were 
computationally expensive so Ward [War92] developed a de- 
vice to acquire the reflectance of a surface sample, as well as 
a mathematical description of the reflectance of anisotropic 
materials. The Ward BRDF uses a Gaussian lobe controlled 
by the standard deviation. For the isotropic version of the 
model this parameter is Om. The model uses the halfway 
vector h and two constants for the combination of the diffuse 
and specular components of the reflection. This BRDF can 
be validated using experimental measurements. The mathe- 


matical expression is: 
2 / (b-n) 
ks exp |- tan ( oe )| 


cos(Wo)cos(w;) ATAR 


friso (Wo, wi) a 


Note that the exponent is zero in the case of the ideal 


reflection. The normalization factor mor’ guarantees the 


correct integration of this function on the hemisphere of 
directions. Because of the Gaussian, the model produces 
quite high values —which should be avoided in a rendering 
system— specially when 0 is around 0.2 at grazing view- 
ing angles [MWM*98]. The behavior of the Ward model at 
grazing angles has been discussed in [NNSK99, GMD10]. 
The use of the normalization term prevents for these numeri- 
cal instabilities, but the given here is not the only. In [Diir06], 
a study on the energy conservation of the Ward BRDF results 
in other possibilities for normalization. 


The anisotropic version of the mathematical formula- 
tion [Wal05] uses an elliptic Gaussian whose axes are 
aligned in a way that is controlled by parameters 01, and Qty. 


frani(Wo, Wi) m 
pie 
k ~ (hen)? 
ATO, Oy y/ cos(Wo) cos(w;) 
where 
(h-x) = sin(Ow, ) cos(Ow, ) + sin(Ow; ) cos(Ow;) 
h 
iy = sin(Ow, ) sin(Ow, ) + sin(Ow, ) sin(Ow;) 
: h 
_ cos(Ow, ) +cos(Ow;) 
(h-n) i 
|h| = [2+2sin(Oy, ) sin(Ow,) (cos(w, ) cos(Ow,; ) 


nie 


sin(ow, ) sin(@w,; )) + 2cos(Ow, ) cos(Ow;)] 


This BRDF is not plausible but is one of the most versa- 
tile reflectance functions as is cheap to evaluate, has a direct 
sampling method and fits well to measured reflectance data. 
Examples of its 3D plot are found in Figure 6 fifth column; 
the fourth column shows examples of our next reflectance 
model. 


Variants have been developed more efficient computa- 
tionally, better suited for fitting measuring reflectance data 
and efficient for Monte-Carlo importance sampling (see sec- 
tion 5). It is the case of the Ward-Diir BRDF model for spec- 
ular anisotropic reflection [GMD10]. 


The Schlick BRDF 


The objective of the Schlick-BRDF [Sch93, Sch94b] is a 
plausible model which is based on the microfacets model, 
the Fresnel equations and the set of parameters also allow 
the distinction between isotropy and anisotropy. Schlick pro- 
poses an alternative way to [SHSL97] of representing he- 
terogeneous surfaces: using a two layer model. The outer 
layer is used for specular reflection and the inner one is 
for diffuse subsurface scattering. Homogeneous surfaces are 
represented by the single layer model with values (C},r, p). 
A heterogeneous surface should use the double layer model 
with values (Cy, 7, p) and (C),r’, p’), where: 


e Q, is the reflectance which depends on A, and gives values 
between 0 and 1. 

e r € [0,1] describes the roughness. When r = 0 it rep- 
resents a mirror; if r = 1 we have the perfect diffuse 
surface. 

e p € [0,1] describes the isotropy. When p = 0 is full 
anisotropy and if p = 1 is full isotropy. 


The meaning of these parameters is the same in both the 
single or double layer model. This implies that the compu- 
tation of S(u) in Cy and S'(u) is the same for C}, only the 
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Ashikhmin BRDF Blinn BRDF 


Poulin-Fournier BRDF 


Schlick BRDF Ward BRDF 


Figure 6: Some examples of 3D plot of BRDF models. From left to right: Ashikminn-Shirley with parameters nu = 1, nv = 75, 


ks = 0.5, kd = 1 (up) and nu = 75, nv = 1, ks 
kd = 0.5, n = 100 (down); Poulin-Fournier d = 2, h 


0.7, kd = 0.7 (down); Blinn ks 


0.3, kd = 0.7, n = 15 (up) and ks = 0.5, 


0.01, n 


10, ks = 0.8, kd = 0.2 (up) and d = 4, h = 0.2, n = 200, 


ks = 0.8, kd = 0.5 (down); Schlick single layer [0.3,0.1,0.2] (up) and double layer {0.05,0.17,0.11][0.12,0.4,0.1] (down). 


parameter value differs. The BRDF model is expressed as: 
fr(Wo, Wi) = R( (n- h), (Wo g h), (wi : n), (Wo i n), on ) 


where the arguments of the reflectance function R are 
substituted by (t,u,v,v’,w) for its calculus. We have 
R(t,u,v,v',w) =S(u,C)) Dit, v,v’,w) in the case of using the 
single layer model and for the double layer it is extended as: 


S(u,C,) D(t, v, v ,w)+ 


R(t,u,v,v ,w) = 7 IN of 1) p/ , 
[I-s (u, C )] S (u,C,,)D (t,v,v ,w) 


Here the specular term S(u,C}) is an approximation of the 
Fresnel term and D(t,v,v’, w) is the directional term that ac- 
counts for emission and re-emission of light at subsurface 
levels. The directional term contemplates the geometry of 
the microfacets G(v,r), and models the anisotropy with a 
dependence on the zenith angle, Z(t,r), and the azimuthal 
angle A(w, p). The directional term is expressed as: 


Gv) Gv’) (1 = Gv) Gv’) 


D(t,v,v',w) = 
(rvv w) 4rvv' 4rvv' 


Z(t) A(w) + 

In addition to providing the previous function, Schlick in- 
troduced in [Sch94a] a computationally effective approxi- 
mation to Fresnel’s equations, which was widely used in 
other models and is present in nowadays applications. 


The Lafortune BRDF 


One of the most multifunctional BRDF models is the one 
introduced by Lafortune et al. [LFTG97]. It was used to fit 
measurements from real surfaces and can be considered as a 
generalization of the Phong BRDF [Pho75] that is reciprocal 
and obeys the energy conservation law. With the reflection 
operator around the normal direction ry = 2(u - n) n — u, 
expressed in a canonical reference system, ry is equal to 
(—xu,—yu,Zu), which is interpreted as a scale factor (-1,- 
1,1) applied to the original vector. When we apply the gene- 
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ralization, a matrix M is used for the reflection operator and 
thus f- œ (w;-Mw.)”. The modified Phong expression is: 


fr(Wo,Wi) = (Wo-tw:)" = (wi-tw,)” 


If we express the diagonal of a matrix M as a vector 0 — 
note that this is a parameter of the model— a retroreflection 
effect is obtained with o set to (1,1,1). By using more than 
one lobe, this BRDF is able to represent the data from ac- 
quired materials. Let N be the number of lobes, which is 
used in the final expression: 


Pa + Lwi (xuX0;; YuYo;; ZuZo;))” 
l=1 


N 
fr(Wo, Wi) = 


4. Function Approximation for Acquired BRDF models 


Suitable analytic models are not always available for a de- 
sired type of surface. Sometimes their formulation does not 
make them useful in a rendering system. To overcome these 
limitations, a number of researchers have proposed and used 
techniques for obtaining BRDF models from measured re- 
flectance data. This strategy allows integration of realism 
in global illumination rendering systems by using data from 
real materials. However raw tabulated reflectance data can- 
not be directly used in rendering systems because it requires 
a considerable amount of memory and processing time, due 
to the high dimensionality of the problem domain. 


A common solution to this problem is to choose a suit- 
able function space, and adapt algorithms to approximate 
measured BRDF functions in that space. In this line, some 
strategies which have been specifically designed include: 


e Basis functions, such as a single Gaussian lobe [War92], 
or summation of cosine lobes [LFTG97]. 
e Spherical harmonics [WAT92, KSS02]. 
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e Weighted sum of separable functions [Fou95, KM99, 
LRRO4]. 

e Zernike polynomials [KvDS96]. 

e Spherical wavelets [LF97,CBP02,CPB03, CPBO06]. 

e Product of functions (factorization) [MAA01, LK02]. 


Acquisition of BRDFs can be carried out using a goniore- 
flectometer [War92,GA97], a device capable of measuring a 
BRDF by moving a light and a sensor through different po- 
sitions in the space of incoming light and outgoing reflection 
(captured by the device), taking measurements every few de- 
grees. After, a denoising preprocessing, a tabulated BRDF 
is obtained. Other researchers have worked towards BRDF 
acquisition methods based on simpler devices, allowing to 
measure material reflectance by using a single digital came- 
ra in combination with curved mirrors [Dan01]. Also, we 
could use computer simulation or emulation of light propa- 
gation in a given type of material, by developing adequate 
models of those materials [CMS87, WAT92, APSOO]. 


After raw data has been obtained (usually a set of re- 
flectance values corresponding to a set of pairs of incoming 
and outgoing directions in the hemisphere), it is necessary 
to filter it to eliminate noise, and then find the closest ap- 
proximation in the given function space. This can be done 
by projection (for function spaces spawned by orthonormal 
basis sets) or using non-linear function fitting techniques for 
other schemes. Usually, the need to save memory implies 
that some compression scheme should be used, finding an 
adequate trade-off between approximation quality and mem- 
ory consumption. 


Fitting is carried out by numerical optimization tech- 
niques which minimize the error between modelled (f’) and 
measured data (f). As stated, measured data is a set of n re- 
flectance values {f|, fo,..., fn} obtained for some set of n 
sample directions pairs, while f’ is a function which can be 
evaluated for those sample directions yielding another set of 
n values { f], f3,..., fa}. The goal of function approxima- 
tion is to find a function f’ which minimizes the distance 
between f’ and f. Two error metrics are commonly used: 


1 a 2 
RMS error o= [nif 


relative RMS error ©; = 


although relative RMS seems more appropriate, we must 
take into account that for this metric it is necessary to ex- 
clude those samples from the summation for which f; = 0, 
even when f # 0. The error metric [NDM05] also could be 
used as the L? metric over the incident and outgoing hemi- 
spheres of directions. This metric has an advantage when 
applied to measurements as it is able to remove the part of 
the domain where data is not reliable. This situation could 
occur with extreme grazing angles (angles larger than 80 de- 
grees) [NDMOS]. 


Function approximation can also be used to approach 


known analytical BRDF models, because this allows to han- 
dle a very wide class of BRDFs by using an uniform ap- 
proach in rendering systems. Fitting exact known analytical 
models is also used to verify validity and performance of 
function approximation strategies. 


Actually measured reflectance data is publicly 
available from various sources, such as the initial 
Ward’s data set [War92], the CUReT (Columbia- 
Utrech) [DGNK99], the Cornell graphics laboratory 
(http://www.graphics.cornell.edu) or the MIT’s MERL 
database [MPBM03]. These data sets differ considerably 
in many aspects, such as the file format or the number of 
samples acquired. In this section we review the literature 
to show the strategies used for BRDF approximation in 
function spaces. 


4.1. Basis Functions 


Specular BRDFs can be approximated by a single Gaus- 
sian lobe, as described in [War92] or by using a summa- 
tion of generalized Phong lobes [LFTG97] (both in sec- 
tion 3.3) or adjusting to other analytical BRDF models such 
as Cook-Torrance, Blinn-Phong or He et al. This type of 
representation consists of fitting some primitive function to 
a physically-based model or to actual reflectance measure- 
ments. The representation must be positive, reciprocal and 
energy conserving. 


When the basis functions are non-linear, a non-linear op- 
timization technique is required to determine the parame- 
ters of the BRDF. In this case we want to minimize the 
mean-square error of the reflectance function, and usually 
the Levenberg-Marquardt technique [GMW81] is applied. 
This is in essence a good option as it needs fewer terms of 
coefficients than other methods, given that a single function 
can capture a complete BRDF over its incident and outgoing 
direction domain, thus it is uniform and compact. In case of 
the generalized Phong lobes [LFTG97], two or three coeffi- 
cients and an exponent are needed. 


This representation possesses many qualitative and quan- 
titative properties. It is able to represent the directional dif- 
fuse reflection: (1) the fading out of the diffuse compo- 
nent for grazing angles, (2) the specular mirror term at gra- 
zing angles, (3) retro-reflection —scattering in the backward 
direction—, and (4) the off-specular reflection —maximum 
directional-diffuse lobe at the grazing direction—. Usually, 
this type of representation omit the Fresnel factor and the fit 
of analytic reflectance data can result in significant error. 


As an alternative, we could use a linear model for approxi- 
mating BRDFs. In [OKBG08], isotropic and anisotropic re- 
flectances can be achieved with a simple and efficient model 
which avoids complicated parameters or non-linear opti- 
mization processes by employing polynomial functions. 


4.2. Separable Decomposition or Factorization 


Factorization techniques represent the BRDF using lower 
dimensional functions, generally as a sum of them, though 
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in the next section we introduce some techniques using the 
product of the terms. Each term is the product of two func- 
tions, one of the incident direction and other of the reflected 
direction. 


Representations based on dimensionality reduction have 
to construct a full sample matrix and perform basis decom- 
position. This process could become significantly costly if 
dense scattered data is used, but is suitable for representing 
measured data with generic analytical models achieving high 
compression rates. By approximating a BRDF table with a 
set of factors we can significantly reduce storage require- 
ments and create a data format more suited to graphics hard- 
ware. The factorization procedure decomposes the BRDF 
function into a sum of terms, generally expressed as: 


N 


fr(wo, wi) = $ ge(Wo)he(wi) (14) 
k=1 


Computation of factors can be done analytically or 
numerically and in the case of real-time rendering they 
are stored in texture maps. There are various ways to 
achieve an analytical decomposition in the sum of pro- 
ducts: from Taylor series expansion to spherical harmo- 
nics [WAT92], Zernike polynomials [KvDS96], K-means 
clustering [LMO1] or as a simple and universal numerical 
approach as singular value decomposition (SVD). 


Given a large enough number of factors, any BRDF can 
be approximated to arbitrary accuracy. However, we need 
a shortened series of factors which should be optimal for 
computational use. For this reason researchers have devised 
many factorization schemes. The numerical algorithm used 
to compute each factor is another source of variation. In the 
following subsections we review some of them. 


4.2.1. Usage of Singular Value Decompostion 


The use of singular value decomposition (SVD) was pre- 
sented in [Fou95, Rus98]. Given a matrix M (of size p x q, 
p > q) this process decomposes it into the product of three 
matrices M = UWV". With the use of SVD each element 
of M is seen as a weighted sum of the rows of U and the 
columns of V”. In the case of the BRDF we take a function 
of four variables and separate them into the sum of the pro- 
ducts of the two two-variable functions: p pairs of (Ow, , Ow, ) 
for the outgoing direction and q pairs of (Ow;,w;) for the 
incoming direction if M; j = f-(Wo, Wi). This decomposition 
obeys reciprocity in [Rus98]. The Phong BRDF and Ward’s 
experimental data was approximated by the authors to be 
used for radiosity calculations. Nonetheless with SVD, the 
existence of positive and negative factors can lead to inade- 
quate results and complicate the parametrization. The pres- 
ence of negative factors also prevents its use for BRDF sam- 
pling purposes, and clampping these values will not solve 
the original problematic. 


4.2.2. Normalized decomposition on graphics hardware 


Kauzt and McCool [KM99] proposed an approximation 
of the BRDF based on textures, expanding the tech- 
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nique of [Fou95]. This is a compressed representation that 
uses naive two-dimensional hardware texture capabilities to 
achieve interactive rates of rendering. In this case the two- 
dimensional functions are held in texture maps. They eva- 
luate two algorithms for finding appropriate functions for g 
and h (terms of Eq. (14)): SVD and normalized decompo- 
sition (ND). The fact that SVD introduces negative terms 
is numerically problematic and is not compatible with the 
strictly positive nature of the reflectance. The normalized 
decomposition is considerably faster than SVD, takes less 
memory and contains only positive factors. 


To improve the BRDF separability, it is recommended 
a change of parametrization to align the main features of 
the BRDF in the matrix representation. The Rusinkiewicz’s 
reparametrization or the elevation/azimuth works well for 
analytical and acquired reflectance data but is incompatible 
with the texture coordinate computations of the hemisphere- 
maps. In real-time rendering, the best solution for hard- 
ware accelerations is to use the Gram-Schmidt half-vector 
parametrization. (See section 2.1 for a description of those 
parametrizations). With this representation it is quite easy to 
evaluate the BRDF at a specific pair of directions, as it re- 
quires few texture lookups, multiplications and additions. 


4.2.3. Product of lower dimensional factors 


The BRDF factorization given by Lawrence et al. [LRR04] 
divides the BRDF function into the product of lower dimen- 
sional factors, stored in a tabular and compact form. Further- 
more, they introduced a novel BRDF sampling procedure. 
The BRDF expression use two distinct factorizations: 


J K 
fr(Wo, Wi) cos(wi) ~ )" Fj(wo) $, Ujx(8w) Vie(dw) (15) 
7 k 


where w is a sampled random vector. The choice of w (the 
parametrization) will depend on the features or complexity 
of the BRDF. In the case of diffuse BRDFs, it is better to 
use w as the local representation of w;. For more specular 
BRDFs, is preferable to use w as h in order to align impor- 
tant properties such as the highlights. 


The initial data matrix Y contains Nw x Nw, samples of 
the BRDF along the outgoing 9 and the outgoing © angles. 
The first factorization, after a reparametrization based on 
the half angle [Rus98], results in the product of 2D factors, 
stored as matrices. Let G be a Nw x J matrix and F be a 
J x Nw, one. The G matrix multiplied by the F matrix is an 
approximation of Y. The second factorization of G, the view 
independent matrix, leads to the product of two 1D matri- 
ces, U and V. Using the non-negative matrix factorization 
(NMF) method [LS00] we always obtain matrices with po- 
sitive values. The resulting L factors, where L = J x K, show 
an intuitive approximation of a specific lobe of the original 
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BRDF. 


Fi (Wo) U; (Ow) Vi (Ow) 


Me 


fr(Wo, Wi) cos(w;) =~ 


II 
Me 


Ti (Wo, Ow, Ow) 


We intended to use these factors as a PDF (Probabilistic Dis- 
tribution Function), thus an upper limit should be calculated 
in order to normalize the factors. This sampling method was 
applied to the Cook-Torrance BRDF model [CT82], Poulin- 
Fournier cylindrical model [PF90], the anisotropic expres- 
sion of Ward [War92], and also to acquired data from the 
public database of Matusik et al. [MPBM03]. Each BRDF 
must be efficiently and independently factorized making this 
method difficult to use for scenes containing many materi- 
als. Some functions will require many factors in order to 
be properly approximated. Thus, the user has to select the 
value of seven parameters manually: J, K, (Nos, x Now)? 


(No, X No, ) and which parametrization to use. 


4.2.4. Directional and spatial variation approximation 


The SVBRDF extends the BRDF to include spatial varia- 
tion in reflectance making this function 6-dimensional. The 
work of Weistroffer et al. [WWHLO7] intends to compactly 
represent data from spatially varying reflectance measure- 
ments, overcoming prior approaches by the fact of using a 
second level of indirection. The decomposition gives a sum 
of products of 2D (the blending weight) and 4D functions 
(the basis BRDF): 


K 
f(x Wowi) = $, Fie(x)Gx(wo, wi) (16) 


The set G of primary basis is defined with respect to a 
secondary linear basis Y: 


K 
fr(x, Wo, Wi) = y F(x) 
k=1 


L 
È TuTi(wo, wi) 
l=1 

where Iy is the weight vector that defines the primary linear 
basis. 


One of the main virtues of this decomposition is the 
choice of the secondary basis. This allows selection of the 
more appropriate representation according to the nature and 
density of the available data. Y could be spherical harmo- 
nics, wavelets (both discussed later in sections 4.4 and 4.5), 
acquired BRDFs or radial basis functions (RBF) [ZREB06]. 
The RBFs are designed to interpolate arbitrary smooth func- 
tions but they need to rely on a larger number of input sam- 
ples. A comparison of RBF with the MERL isotropic BRDF 
dataset [MPBM03] concludes that the more restrictive basis 
of a dataset (the acquired one) is preferable, given that it is 
meant to represent reflectances. 


4.2.5. Finding the required complexity 


Mahajan et al. [MTRO8] gave a theoretical and empirical 
analysis of the BRDF factorization. They answer an impor- 


tant question: how many terms do we need to represent a 
BRDF model accurately with a given parametrization? In 
some cases, as the Lawrence et al. factorization [LRR04], 
it is not just a single question. The user needs to decide 
which parametrization to use, the size of the initial data and 
of course, the number of factors. Mahajan et al. analyze the 
Phong BRDF starting from its spherical harmonic expan- 
sion, concluding the following statement: the numbers of 
terms needed are linear in the Phong exponent, and quadratic 
in the frequency content along the reflected or half-angle di- 
rection. With experimentation they found that this affirma- 
tion holds quite generally for many BRDFs from the MERL 
isotropic database [MPBM03]. 


4.3. Product of terms 


The Homomorphic Factorization (HF) [MAA01] decom- 
poses arbitrary BRDFs into the product of two or more 
factors of lower dimensionality and positive values. This 
method does not require tabulating the data into a matrix of 
the same resolution as the final representation. It deals with 
a sparse data matrix without resampling or interpolation. In 
this case, the decomposition is not a sum of terms but a pro- 
duct of positive factors of two dimensions. Those terms can 
be stored in texture maps after projection with function 7. 


J 
fr(Wo, Wi) © [4 (aj (wo, wi)) 
j=l 
Then the authors use a logarithmic transformation of both 
sides of the equation. The resulting system is a linear data- 
fitting problem. An interesting aspect of this logarithmic 
transformation is the fact that minimizing the RMS error 
in logarithmic space is the same as minimizing the relative 
RMS error in the original space. This improves the solution 
as the perception of the human visual system is relative to 
the ratio of intensities and not to absolute quantities. 


J 
log fr(Wo, Wi) x y log tj(1j(Wo, Wi)) 
j=l 
The resolution of this system has to be reconverted to the 
original problem by exponentiating. Special care has to be 
taken with the logarithmic of zero-value reflectance data. 
Fixing J, we have a three-factor approximation of the BRDF: 


fr(Wo, Wi) © F (wo) G(h) F(w;) (17) 


log fr(Wo, wi) © log F(wo) + log G(h) + log F(w;) 


The numerical technique employed to find F and G is a li- 
near system of constraints relating each sample of a BRDF 
to the corresponding texel locations in F and G. Additional 
equations are included in the linear system to set smooth- 
ing and weighting constraints to the solution. The homo- 
morphic factorization expression Eq. (17) can be used with 
isotropic and anisotropic BRDFs. This representation is si- 
milar to SVD but avoids negative numbers and its flexibility 
allows to use any parametrization. Unfortunately McCool et 
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al. method does not allow for proof of optimality. A study 
and improvement of HF is given by [Col02] whose technique 
automatically selects the projection functions. 


Another important BRDF approximation approach was 
presented based on the HF. In [LK02] the factorization con- 
siders the product of the BRDF and the lighting. Using a 
product function instead of a function, is in fact more com- 
plex, given that lighting uses global parametrization (and is 
defined over the sphere of directions) whereas the BRDF is 
defined with local parameters. This complexity is resolved 
by the use of isotropic BRDFs in arbitrary lighting condi- 
tions and a preprocessing step. With an isotopic BRDF not 
only do we have a reduction in the dimensionality (now 3D), 
but more importantly, the local coordinate frame n,t,t x n 
does not depend on the tangent vector. 


The HF takes the logarithmic of both sides, yielding 
log Ly (Wo) ~ log tı (n) + log t2 (rw, ). It is solved as a linear 
equation system with the Quasi-minimal residual (QMR) al- 
gorithm. The final solution comes from the exponentiation. 
As in [KM99] the approximation is based on textures, and 
the number of factors is fixed to J = 2, so only two textures 
are implemented. Using a mapping 7, texture coordinates are 
converted to directions and vice versa. 


4.4. Spherical harmonics 


Spherical harmonics (SH) are the Fourier basis functions ap- 
plied on the sphere, which is an upper set of the directional 
hemispherical domain of the BRDF. The use of a global re- 
construction scheme such as the SH basis functions has been 
used in the literature of Computer Graphics. Some notewor- 
thy papers are listed here: 


e Cabral et al. [CMS87] have the idea of using the SH to 
derive isotropic BRDFs. They also give an expression of 
the lighting integral which was a simplified dot product. 
To reach that point, they make some simplifications: (1) it 
is assumed constant view direction, and (2) the SH coef- 
ficients are derived from a height field geometry, neither 
from an analytic evaluation of a BRDF model, nor from 
measured reflectance data. 

e Westin et al. [WAT92] use spherical harmonics in the pre- 
computational inference of a BRDF from geometric mo- 
dels. The view and lighting dependence on the BRDF is 
encoded in a matrix that stores the SH basis coefficients. 
In this case the evaluation of high-order basis functions is 
undertaken during rendering, and that implies low com- 
putational cost. 

e Kautz et al. [KSS02] represent the four dimensional space 
of an arbitrary BRDF (any isotropic or anisotropic re- 
flectance model) in a 2D table of spherical harmonics co- 
efficients. Each entry represents the product of the inci- 
dent light times the BRDF product function (the BRDF 
multiplied by the cosine of the normal with the incident 


Montes & Ureña. University of Granada 


direction). This representation is intended to account for 
low-frequency lighting and arbitrary BRDFs in real-time. 
This approach needs to compute the BRDF and the ligh- 
ting coefficients for the Spherical Harmonic basis. The 
first one is precomputed by numerical integration. They 
first parametrize the BRDF product function by local view 
direction, which is represented by SH basis. 


N 
r(Wo, Wi) cos(w;) = L <i(wo)vi(wi) 


ci(Wo) = yi(wi) fr(Wo, Wi) cos(wi) max(0, zv) do(wi) 


The above integral —called the SH-projection— is also 
applied to the incoming lighting. By representing the 
lighting function as well as the SH basis, the lighting inte- 
gral is reduced to a simple sum of dot products [CMS87]. 
To perform integration, we need to rotate the lighting into 
the BRDF’s local frame with a single computation in SH 
space. The rotated coefficients of c; are computed via li- 
near transformation. Moreover, for low-frequency ligh- 
ting, few components (N = 25) are required and the com- 
putation of this integral and its projection into SH basis 
can be performed on-the-fly. If both lighting and mate- 
rial contain higher frequencies, then the number of coeffi- 
cients (understood as resolution or quality of the approxi- 
mation) should be incremented to achieve good results. 


4.5. Wavelets 


If we have decided to represent the BRDF in tabular form, it 
can be helpful to build on the model of wavelets given that 
this approach is widely employed in other fields of science 
for data compression (in particular 2D images). The wavelet 
transform WV is a linear operator that can encode a general 
function f in a process named the analysis of the function f: 


W:(f:4— B) > T(P) 


T (af +Bg) = aT (f) + BT (8) 


The analysis consists of projecting the function at different 
resolutions or scales / € N. The mother basis function Y; = 
yi’ : A — B is the wavelet itself. Commonly the Haar basis 
is used as mother function. The scaling functions ®;}; = 
of 41 :A — B are used to derive the wavelets at resolution /. 
The wavelet transform is a recursive decomposition of the 
function f at different levels or resolutions from the coarse 
level / = 0 to the finest l = n. 
n 
f= Yat Y arw 
k 1=0 m 

where ak and dj" are the projection coefficients. The dis- 
crete wavelet transformation produces the same number of 
coefficients as the original samples, but those whose value is 
zero or below a threshold are discarded and not stored, thus 


achieving compression. The synthesis or reconstruction of f 
uses a similar recursive procedure. 


The more outstanding qualities of this representation are: 
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e generality: it can represent any function, 

e compression: it can achieve important compression rates 
most notably with large datasets, 

e trade-off between accuracy/space easily controlled, 

e multiresolution: reconstruction at different levels of accu- 
racy, 

e reconstruction from the piece-wise function is effectively 
computed in O(logN) where N is the number of BRDF 
samples, 

e simple error metrics, 

e denoising: the wavelet smoothly interpolates the original 
data, 

e it preserves high and low frequency features during 
the decomposition step, and therefore typical reflectance 
peaks of a glossy BRDF are well preserved, 

e there are sampling procedures to use the approximated 
BRDF in a Monte-Carlo algorithm for Global Illumina- 
tion. 


Some of the most representative works related to BRDF 
compression using wavelets are commented upon in the fol- 
lowing subsections. 


4.5.1. Spherical Wavelet Decomposition 


The more common wavelet functions are on R but there 
are also Spherical Wavelets [SS95]. This basis constructs 
a regular subdivision of an octahedron or a tetrahedron 
(for S? or Q spaces) which is independent of the spheri- 
cal parametrization of the incoming and outgoing directions. 
After the subdivision of the domain in triangles, the BRDF 
is approximated by a piece-wise constant function on the tri- 
angle at a specified level of subdivision. The value of the 
function —a single value Fr (Wo, Wi) or the whole spectrum 
response f;-(Wo,Ww;,A) in [CBP02]— is stored in the trian- 
gle in the direction of the triangle centre. The scaling and 
wavelet functions are an area extension of the Haar basis 
functions, and they are known as Bio-Haar basis [SS95]. 


4.5.2. Multidimensional Wavelet Decomposition 


Lalonde and Fournier [LF97] presented an approximated 
BRDF model based on the projection of the BRDF onto 
a 4D wavelet basis. This dimensionality comes from the 
product of four 1D basis that are computed following the 
standard decomposition —the non-standard approach would 
need to use the product of multidimensional basis— after a 
reparametrization of Q x Q domain into R*. More specifi- 
cally (8,) € Q is mapped to (x,y) € R?. 


The BRDF is decomposed as follows: 


fr(Wo, Wi) = Fir (Kwo s Ywa: Kwi, Yw) 


= EEEL Bel, )B (ns); (om BCH 


E hj 

The general bases are defined as: 
(x) if n >0 with l,m >0 
Baa) = 4 ryt À l 
2/°Y(x—m) if n=2 +m 


where ¥ is the mother wavelet, ® is the smoothing function 
and / is the level in the wavelet subspace. 


Favorably, this representation includes the cosine term for 
the projected area. One minor drawback of this represen- 
tation comes from the parametrization, which filled corner 
zones with zeros and they are stored even though they give 
no information. If there are so many wavelet coefficients as 
original data points, we have no compression at all. For this 
reason, the selection of a convenient data structure is im- 
portant. Lalonde and Fournier recommend the use of a hash 
table or a wavelet coefficient tree that stores only the coeffi- 
cients indexed by the wavelet basis selector, however intro- 
ducing a pointer set cost. 


To reconstruct the BRDF at a given set of directions d = 
(Kw, , Yw, » Kw;, Yw; ), each coefficient for which Bn is non-zero 
needs to be consulted. Resulting time complexity is o(w*i*) 
where w is the wavelet width. 


4.5.3. Generic Wavelet Decomposition 


This is a wavelet representation for a reflectance function de- 
pendent on the wavelength of incoming light —also known 
as spectral BRDF— and thus with five dimensional com- 
plexity. Furthermore this representation is also suitable for 
path tracing. It appears in [CPB03, CBP, CPB06] and is 
based on the separation of the BRDF into two components: 
a spectral and a geometrical component. The function fr is 
then Q x Q x R > R. Fixing wo, we have f;(w;,A) : Q x 
R —> R which is equivalent to f-(w;) :Q—> (f-(A): RR). 
This transformation is key for the generic wavelet transform 
because now we have a spherical function f,(w;) whose ele- 
ments are one-dimensional real f(A) functions. The authors 
apply multiple simple transformations. Firstly a spherical 
wavelet changes the spherical signal and gives as a result 
the wavelet and scaling coefficient vectors. Secondly a one- 
dimensional transform is applied to vector data resulting in 
the final coefficients. The second transformation is more ef- 
fective and compresses more than the spherical domain be- 
cause spectral values tend to be smoother. The combination 
of the spherical transformation and the spectral one provide 
the full transformation of the spectral reflectance distribution 
function, obtaining a compression rate of more than 90%. 


As we have already pointed out, coefficients with zero or 
a value below a threshold are ignored, leaving a sparse set 
of data which has to be stored efficiently. In this case the 
recommended data structure [CPB03] is the sparse array, a 
double-linked list of strips. The strip node stores a bit ar- 
ray of valid/invalid positions together with the array data. 
An optimal number of strips is crucial for the storage sa- 
ving of this structure. In [CPB06] a two-threshold scheme 
is recommended, one for the spherical transformation and 
other for the spectral transformed coefficients. The fact that 
the BRDF diffuse component yields more compressed re- 
sults even by using different thresholds for each of the BRDF 
components, makes them to develop a formula for adaptively 
computing the threshold. This local value ensures a relative 
uniform compression. 
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In [CPB06, CBP07] a general Generic Wavelet Transform 
(GWT) framework is intended to be used with acquired 
spectral reflectance data in a real-time per-pixel rendering 
system. The coefficients of the 4D wavelet transformation 
are stored in 3D texture maps and their processing is per- 
formed on the GPU using a fragment shader. The full 4D 
decomposition is done in two steps. First they transform the 
sphere domain S? in R? with a given mapping. Then, the 
spherical wavelet transform [SS95] is applied twice: the 2D 
wavelet transformation for the incoming and the outgoing 
hemispheres of directions. Directions correspond with tri- 
angles and wavelets, and each triangle data is stored on a 
high-resolution texture. The number and textures that the 
fragment shader uses depend on a resolution parameter that 
is chosen automatically for the BRDF reconstruction. This 
level-of-detail reconstruction allows for optimal antialiasing 
of surfaces at a distance from the observer. 


5. Sampling the BRDF 


The way of distributing random samples is known as the 
probability distribution function (PDF). In Monte Carlo al- 
gorithms for Global Illumination the generation of random 
sample directions (to estimate the reflected radiance) is a fre- 
quent task and should be performed as quickly as possible. 
In order to render realistic images, we have to use realistic 
material models, although few of them provide a method for 
sampling with the desired PDF and properties (that is, pro- 
portional to the BRDF and fast). Some BRDFs [TS67,CT82] 
are impractical in rendering systems for this reason, tough 
some approximations could be used as the ones presented in 
the following sections. 


A good definition of the PDF is the key to the effective- 
ness of the final MC estimator. It is preferable to generate 
more samples where the function has higher values —that 
is what importance sampling (IS) does— because with uni- 
form sampling we will have under-sample regions in which 
high values of the function occur, leading to visible error 
in the image. It is also important to ensure that a region 
with non-null BRDF values would have a non-null proba- 
bility to be sampled. Importance sampling gives better so- 
lutions when the PDF is chosen closer to the integrand but 
it is not easy to find for a wide-range of BRDF models. A 
possible solution to use the BRDF for IS is to have a math- 
ematical framework that forces plausible BRDFs. Edwards 
et al [EBJ*06] create an empirical BRDF from the trans- 
formation of halfway vectors into different domains assum- 
ing its integral can be a PDF. This constraint enforces en- 
ergy conservation. Thus the unit radius disk is a good do- 
main to specify BRDF properties such as specular and retro- 
reflection. 


5.1. Monte-Carlo based BRDF sampling 


If we consider any algorithm which generates random direc- 
tion vectors (or samples) in the hemisphere Q, the distribu- 
tion of those samples is determined by a probability mea- 
sure, which obviously depends directly on that algorithm. 
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We will symbolize probability as P. For any region R C Q, 
the real value P(R) is the probability for the random vector 
to be in region R. 


Measure P usually depends on a vector v in the hemi- 
sphere. For clarity, we will express this dependence expli- 
citly in the notation by writing Py for the particular probabil- 
ity measure associated to any vector v € Q. Usually Py(Q) = 
1, although for some algorithms, samples are produced in the 
hemisphere under the surface. For those cases, measure P is 
defined on the whole sphere S°, with P(S?) = 1. 


The probability distribution function (PDF) is the density 
of Pu with respect to solid angle measure in the hemisphere, 
that is for any u,v € Q, the real value py(u) is defined as 
dPy(u) /do(u). When a PDF p is used to sample a BRDF f,, 
it must hold this condition: 


fr(u,v) > 0 = > pu(v) > 0 and py(u) > 0 


When using a PDF py for sampling a BRDF f,, we can 
obtain an estimation of the reflected radiance at u by using a 
weighted average of the estimated radiance values for each 
sample direction {v),...,Vn}: 


1 Š , fr(u, v;i) cos(v;) 

L,(u) F 2 w;L;(v;) where wi mo 

The values L, (v;) can be approximated by using the same 
method recursively, as it is done in raw path-tracing [Kaj86], 
or by using a different method, as in the method proposed 
by Ward et al. [WRC88] (where pre-computed radiance es- 
timates are used) or in final gathering for Photon-mapping 
[Jen01] (where radiance estimation for secondary rays are 
obtained from the photon-map). 


In the following sections we review particular sampling 
strategies for particular BRDF models. In all of the cases, 
the problem is to find a particular map M with the property 
described above, and which can be efficiently evaluated. 


5.2. Direct Sampling 


Usually the PDF is chosen to mimic a specific model, and 
this situation is known as direct sampling. 


Ideal Diffuse Case. A perfect diffuse BRDF reflects ener- 
gy equally distributed in Q. This simple BRDF (f, = 
kq/™, where 0 < kg < 1) has an associated PDF pu(v) = 
cos(v)/1 which is exactly proportional to the BRDF times 
the cosine term and also independent of u. To get a ran- 
dom vector v distributed according the PDF, compute: 


(8y,0v) = (arccos( /&1), 2262 ) 


where €, and €) are two independent uniformly dis- 
tributed random variables with values in [0, 1). From this 
point we will use & notation with this meaning. 

Ideal Specular Case. The PDF is a Dirac-delta function 
Pu(V) = Ou(v), and thus a sample distributed according 
to this PDF is ry with probability 1, and any other values 
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have probability 0. As a consequence, generating random 
directions consists simply of producing ry for a given u. 

Lobe Distribution Sampling. BRDFs based on cosine- 
lobes [Pho75, Bli77, Lew94, LFTG97, War92] have an as- 
sociated algorithm for sampling. The single-lobe BRDF 
is defined as: 


fr(u,v) = C(n) (v: ru)” 


where n > 0 is a parameter, and C(n) is a normalization 
factor which usually depends on n. This ensures these 
BRDFs have an energy conservation property. 

For this BRDF class, the directional-hemispherical reflec- 
tivity can be written as: 


Pn(u) = C(n)N2(tu,n) 
where, for any vector a, M2 (a,n) is defined as: 


def 


Nolan) È f (v-a)" (v-n)do(v) 


which is called a double axis moment. Arvo [Arv95b] pro- 
vided an analytical expression for M2. It is easy to show 
that M2 (a,n) < N2(n,n), as Na (n,n) = 27/(n + 2) (thus 
it can be easily evaluated). A usual option for C(n) is 
1/N2(n,n) ensuring energy conservation and yielding a 
simple BRDF model which can be evaluated very fast: 


+2 
fluv) = A Vra)” 
A related normalized PDF can be defined as: 
_ 1 n 
Pu(v) = i(Tu,n) (v ru) 


where N ensures normalization and is defined as: 


Mian = 1. (v-a)"do(v) 


This is called a single axis moment and, as for No, it is 
possible to give analytical expressions for Nj. In order to 
obtain samples distributed according to this PDF, we first 
obtain a random vector v whose spherical coordinates are: 


(0,6) = (arccos ( P) ant) 


ry = (sin(8) cos(o),sin(®) sin(), cos(®)) is computed, 
and after this, we rotate ry by using R, where R is the ro- 
tation which converts n into ru, that is ru = R(n). The re- 
sulting sample is then s = R(ry). Note that this yields unit 
vectors outside Q (with a probability smaller than 1/2), 
and in these cases we must reject those samples and try 
to produce a new one until a valid one is obtained. The 
exponent n is a parameter of this sampling method. 

Half-angle based The reflectance models [Bli77, War92, 
KSKO1] calculate vector h as (Oh, 0n) = ( "4/61, 2762). 
Reflected vector v is returned with a probability of 

n+ 1 cos”(Op) 


2n 4(v-h) 


Pulv) = 


The above formulation includes the correct proportion be- 
tween the measures: 


pulv) = palv) g | a8) 
do(h) _ sin(On) dOh dOh _ sin(On) dOh don 
do(u) sin(Oy) dOy doy sin(20,) 2409p dOn 
_ sin(9n) = 1 
4 cos(6,) sin(On) 4 cos(8p) 
1 1 
~ 4(v-h)  4(u-h) 


5.3. General Sampling 


General sampling methods deal with arbitrary reflectance 
models and measured data. This will permit the use of com- 
plex reflectance distribution functions in Global Illumination 
applications and in real-time rendering. Many papers have 
been published for this topic mostly based on sampling by 
inversion procedure, or on the rejection sampling technique. 


e Lalonde and Fournier [LF97] use numerical inversion 
of the BRDF for sampling purposes using the wavelet 
representation. Samples are generated directly from the 
wavelet data structure. 

e The Homomorphic Factorization [MAA01] decomposes 
arbitrary BRDFs or measured data into the product of two 
or more factors which are stored as texture maps. Fixing 
the direction Wo, either G(h) or F(w;) (see section 4.3) 
can be used as a PDF to generate randomly outgoing sam- 
pled directions. If samples are generated according to the 
half-vector some corrections should be done to the final 
distribution, as we have already pointed out in Eq. (18). 

e The MERL BRDF database ' is an approximate represen- 
tation of a hundred isotropic surfaces, with a resolution 
of 90 x 180 x 90. Matusik et al. [MPBM03] developed a 
generic sampling algorithm for these measured data. The 
procedure computes a 2D cumulative distribution function 
(CDF) over incident directions according to the spheri- 
cal parametrization of the hemisphere. This requires dense 
sampling along all variables (even more storage than the 
BRDF itself). Given that a single isotropic BRDF is 33 
MB, this sampling scheme is not compact. 

e The Spherical Wavelet BRDF representation has the ad- 
vantage of a wavelet sampling technique [CPB03]. An 
alternative tree structure independent of the represen- 
tation stores the BRDF fÀ = J fr(wo,wi,A)dX as the 
monochromatic version with the cosine term included. 
To randomly sample a direction s, first they choose a piece 
of the representation, a triangle, using its CDF. Thus, a 
triangle projected onto the hemisphere on level j, T} , is 
then uniformly sampled. The fact that each triangle value 
at level j is the average of its children values at level j+ 1 


T Available on line at: <http://www.merl.com/brdf/> 


Montes & Ureña. University of Granada 


Montes & Ureña / An overview of BRDF models 21 


enables a recursive search procedure making the algo- 
rithm’s complexity dependent on the number of triangles. 
The conditional probability is computed as the probabil- 
ity of selecting a triangle times the probability of the uni- 
formly sampled value. To reduce the bias in their sampling 
method [CBP07] the selection and descent of a child tri- 
angle use random permutation of the children instead of 
a sequential search. The complexity of this algorithm is 
logarithmic in the number of spherical triangles. 

e Wavelet Importance Sampling (WIS) [CJAMJO5] gene- 
rates values that are proportional to the product of two N- 
dimensional functions, the BRDF and the lighting, given 
their representation using Haar-based wavelets. The pro- 
duct is a new wavelet defined by the product coefficients 
of both wavelet representations. The algorithm starts from 
a coarse level of representation towards the fine level of 
the wavelet tree, warping points according to f,- Li by fol- 
lowing a probability defined by the scalar coefficients of 
the child nodes. It is not an exact technique for sampling 
because the representation of both functions, reflectance 
and lighting, comes from sparse levels of detail. 

e The factorization of analytical and measured reflectance 
models [LRR04] includes a sampling procedure. The non- 
negative matrix factorization (NMF), ensures that the fac- 
tors are always positive, and facilitates the use of norma- 
lized factors as a CDF for numerical inversion sampling. 
First the algorithm selects the lobe / that contributes most 
energy for the current view u, according to the F ma- 
trix (see section 4.2.3). The CDF for this sampling is re- 
computed when the outgoing direction changes. Next, the 
hemisphere is sampled according to the selected lobe / by 
sequential generation of the elevation and the azimuthal 
angles using the precomputed CDF for U; and V;. The 
density function for the generated outgoing direction is: 


do(w) 


1 & (u, 0w,0w) 
do(v) £ 


| ~ a(n) yh Fj) 


where w is the parametrization vector, L is the number of 
lobes, and u, v usually stand for Wo and w; respectively. 

e The Cascade CDF method [LRROS] is an adaptive tech- 
nique oriented to sampling acquired reflectance data. It is 
based on sampling by inversion of the CDF, and instead 
of uniformly distributing the samples, it uses a second 
and equivalent distribution which is compact. The method 
starts with a N-dimensional PDF and divides it into the 
product of a 1D marginal distribution ø and a set of 1D 
conditional distributions. 


pal) = aw) | 


J: 
1 T P(x y) dx 
Xi—Xi—1 Ja P’) 

Compression is carried out using the Douglas-Peucker 
greedy algorithm [Ros97] which approximates a curve 
(in this case the CDF) employing an optimal number of 
segments. Applying the aforementioned to BRDF sam- 


plx) = 
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pling the authors recommend starting from an initial dense 
set of data. The required processing time and memory to 
compute the adaptive representation may vary with each 
BRDF. We applied this technique to the MERL BRDF 
database to test the compression of data. Some of the re- 
sults are shown in table 2. With a resolution (Ou x ou x 
Oh X Op) of 32 x 16 x 256 x 32, each uniformly-sampled 
CDF occupies 33 MB. We can see that average compres- 
sion rate is 98%, mostly in non-specular reflectance data. 

e In [MULGO08] an importance sampling algorithm for arbi- 
trary BRDFs and measured data is presented. This method 
has the benefit of optimal rejection sampling, in the sense 
of bounded trials. It has few and intuitive parameters that 
could be set as constants. A scene with many reflectance 
models can be rendered using the same algorithm’s con- 
figuration for each BRDF instance, avoiding user guid- 
ance in terms of adjusting parameters for each BRDF. 
The core data of this approach is a set of hierarchical 
quadtree structures: one for each (f;, Wo) pair. Each pre- 
computed quadtree subdivides the disc domain [—1,+1]°, 
adaptively in a recursive way. Samples on the disc will 
follow a distribution which is exactly proportional to the 
BRDF times the cosine term. The subdivision criteria is 
a parameter of the method which ensures that rejection 
sampling on leaf nodes can be done with an a priori 
bounded number of average trials. This algorithm gives 
favourable results when used with a set of common and 
complex BRDF models and also with MERL isotropic 
BRDF database, even in the same scene. This PDF can 
be integrated with other techniques such as Resampling 
Importance Sampling (RIS) [TCEOS]. 


To compare direct and general sampling we have used 
the Dragon model from Stanford University’. The re- 
flectance function used in this scene corresponds to Oren- 
Nayar [ON94] with a roughness value of 0.83, and an in- 
stance of Strauss BRDF [Str90] mostly smooth for floor 
and wall, respectively. The dragon itself has a Lafortune 
BRDF [LW94] with exponent n = 46. The rendering is in 
figure 7. With this mixture of BRDFs, we can compare the 
Disc sampling method with uniform sampling, cosine lobe 
in Q, and Lawrence et al. [LRR04] factored representation. 
For the cosine lobe and the factored BRDF we adjusted the 
parametrization manually to fit the shape of each BRDF in- 
stance. Note that the cosine lobe technique is suited for lobe- 
based BRDFs and no others. From left to right, the sampling 
time are: 9.749, 29.06, 24.539 and 42.632 seconds respec- 
tively. With only 100 samples, Disc PDF [MULG08] gives 
less noise than the others, with no need for manual adjusting. 


6. Conclusions 


Our main objective is to divulge common understanding of 
a wide range reflectance models. For this reason, we have 


f The Stanford 3D Scanning Repository is on-line at 
<http://graphics.stanford.edu/data/3Dscanrep/> 


22 Montes & Ureña / An overview of BRDF models 


Figure 7: From left to right, images corresponding to the Uniform PDF, adjusted cosine-lobe strategy in Q, the Factored 
representation of the BRDF and the Disc PDF. With only 100 samples, Disc PDF algorithm performs better, as the scene has 


three different instances of BRDFs. 


MERL BRDFs Time (sec) Mem. (KB) Reduc. 
alum-bronze 572.13 998.15 97.05% 
alumina-oxide 577.4 174.75 99.48% 
aluminium 568.31 1156.5 96.58% 
aventurnine 562.68 331.65 99.02% 
beige-fabric 559.53 187.66 99.44% 
black-fabric 579.51 265.13 99.22% 
black-obsidian 580.54 802.2 97.63% 
black-oxidized-steel | 575.38 761.45 97.75% 
black-phenolic 577.38 901.96 97.33% 
black-soft-plastic 562.24 799.8 97.63% 
blue-acrylic 576.79 382.53 98.87% 
blue-fabric 578.26 351.58 98.96% 
blue-metallic-paint2 | 572.29 1038.0 96.93% 
blue-metallic-paint 575.35 1110.05 96.72% 
nickel 567.49 1023.62 96.97% 
red-plastic 586.5 275.83 99.18% 
teflon 561.78 184.51 99.45% 
violet-acrylic 578.69 826.69 97.55% 
white-marble 565.71 247.79 99.27% 
yellow-paint 566.06 143.47 99.58% 


Table 2: Requirements of Cascade CDF method when the original 
CDF is compressed and adaptively sampled. 


presented the fundamentals of the reflectance function, its 
characteristics and the environment in which it is used. With 
the intention of be a complement of information on how to 
choose, use or even implement a specific model we have 
overview a set of papers that cover a wide range of the full 
bibliography. For each BRDF we briefly expose its formula- 
tion and applicability, and we have covered issues of interest 
in representations and sampling methods. 


We also have reviewed important contributions about the 
possible representations of the BRDF and how to sample ac- 
cording to it, giving a survey of techniques to Importance 
Sample the BRDF function in a Monte-Carlo rendering sys- 
tem. 
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